@ChrisRackauckas TaylorDiff.jl works like a charm, although at the moment it does not work with Double64 numbers. I don't completely understand the TaylorDiff.jl API, am I using this right? What does the second to last argument mean?

https://github.com/SouthEndMusic/Integrals.jl/blob/b60c6a467ccb8fce4235c5a4422f89651e73c1d2/ext/IntegralsGaussTuranExt.jl#L91-L99

Edit: l turns out to be the direction of the derivative.

Edit: Using Double64 works when DoubleFloats is imported before TaylorDiff.

Very nice! Just a thought, but have you considered opening a new package for your code to compute the Gauss Turán rules? Packages like GeneralizedGauss.jl, QuadGK.jl, and FastGaussQuadrature.jl are all standalone as they solve a very specific problem, and Integrals.jl isn't as much of a place to do something novel as it is a place to provide a consistent interface to other packages as well as AD rules. This also simplifies some aspects of dependency management and writing tests, and you could still implement the Gauss Turán algorithm as a package extension also depending on TaylorDiff.jl, as you have done here. I'd be happy to advise.

@lxvm That sounds like a good idea! Indeed to evaluate an existing Gauss-Turán quadrature rule you only need TaylorDiff so that nicely fits in an IntegralsTaylorDiffExt.jl. Say I make a package GaussTuran.jl, then I have some questions:

Would GaussTuran.jl need to be a weak dependency of Integrals.jl at all? Or would IntegralsTaylorDiffExt.jl just tabulate computed GaussTuran quadratures and an interface to evaluate them? You at least want to prevent that people need the dependencies of GausTuran.jl (Optim, PreallocationTools) just to evaluate computed Gauss-Turán quadrature rules.
Or would you say that the main purpose of GaussTuran.jl is to tabulate the quadrature rules, and computing your own rules should be an extension of that package?

I think it makes sense for GaussTuran.jl to tabulate quadrature rules, and a good example for how to do this is QuadGK.cachedrule. If a rule is not tabulated, then it could be calculated on the fly. It's up to you whether PreallocationTools.jl and Optim.jl would be (weak) dependencies of GaussTuran.jl, since the tradeoff is having fewer dependencies for tabulated rules versus the convenience of being able to compute the rules without loading additional packages.

As for how TaylorDiff.jl plays into this, there are two options:

• If GaussTuran.jl just computes rules, then an IntegralsGaussTuranTaylorDiffExt would be needed
• TaylorDiff.jl could also be a (weak) dependency of GaussTuran.jl if that package can implement its own interface for computing the integrals, and a IntegralsGaussTuranExt would suffice

I think the latter would make sense so that GaussTuran.jl is a standalone package, similar to QuadGK.jl, and the extension in Integrals.jl would call that interface. The advantage is then that other people can use your package as-is without Integrals.jl, even though you develop the package with Integrals.jl in mind.

@lxvm Any thoughts on which rules are interesting/relevant to tabulate? I have the not very well thought out impression that it is most interesting for sets of functions where the derivatives are not in the span of the set.

You could start with the n=5,s=1 and n=5,s=3 polynomial rules from Ref. 1 in #245. A good example of how to cache rules is QuadGK.cachedrule.

I haven't seen these quadratures in the wild and I think they aren't common because calculating the rule is expensive, calculating higher-order derivatives is non-trivial, and for a given order of accuracy you can trade higher order derivatives for more integration points. With packages like FastGaussQuadrature.jl, getting high-order Gauss-Legendre rules looks very easy and unless there are clear reasons why fewer evaluation points with derivative information is needed, I suppose most people choose a higher-order Gauss-Legendre rule.

I'm starting to get my ducks in a row here: https://southendmusic.github.io/GaussTuranQuadrature.jl/dev/. Let me know what you think 😃

As mentioned by Chris, using Gauss-Turán rules is probably beneficial if you can compute the derivatives of your integrand very cheaply, possibly leading to fewer integrand calls to obtain the same accuracy as with a 'classical' rule with more nodes. See for instance my example: https://southendmusic.github.io/GaussTuranQuadrature.jl/dev/evaluate_rule/.